function out = inflVarianceRatio(outEstim,xSim)

model = outEstim.model;
setupFilter = outEstim.setupFilter;

% Solve for expected inflation up to max maturity+1
vectorMom3   = zeros(1,model.ne);
posInfl= find(strcmp(model.labely,'$\pi_t$'));
maxMat = max(setupFilter.maturities);
% The arrays contain expectations from 1,2,..., Max_maturity into the future
[E_paix,E_paixx,E_paiss,E_paixxx,E_paissx] = CondMoments_Mom3_log(model.gx,model.gxx,model.gss,...
        model.gxxx,model.gssx,model.gsss,...
        model.hx,model.hxx,model.hss,...
        model.hxxx,model.hssx,model.hsss,model.eta,...
        vectorMom3,posInfl,maxMat+1,setupFilter.appOrder);
E_pai0 = repmat(model.g0(posInfl,1),maxMat+1,1);
%[E_pai0,E_paix,E_paixx,E_paixxx]= projectionStepExpectedControl(var0,varx,varxx,varxxx,...
%   maxMat+1,model,setupFilter.setupProj);

% Solve for expected policy rate up to max maturity
posr= find(strcmp(model.labely,'$r_t$'));
[E_rx,E_rxx,E_rss,E_rxxx,E_rssx]=CondMoments_Mom3_log(model.gx,model.gxx,model.gss,...
        model.gxxx,model.gssx,model.gsss,...
        model.hx,model.hxx,model.hss,...
        model.hxxx,model.hssx,model.hsss,model.eta,...
        vectorMom3,posr,maxMat,setupFilter.appOrder);
E_r0  = repmat(model.g0(posr,1),maxMat,1);
% Adding r_cu to the solution
E_r0  = [model.g0(posr,1)      ;E_r0];
E_rx  = [model.gx(posr,:)      ;E_rx];
E_rxx = [model.gxx(posr,:,:)   ;E_rxx];
E_rss = [model.gss(posr,1)     ;E_rss];
E_rxxx= [model.gxxx(posr,:,:,:);E_rxxx];
E_rssx= [model.gssx(posr,:)    ;E_rssx];

% One-step ahead expectations of nominal bond yields
nx       = model.nx;
E1_by0   = zeros(maxMat,1);
E1_byx   = zeros(maxMat,nx);
E1_byxx  = zeros(maxMat,nx,nx);
E1_byss  = zeros(maxMat,1);
E1_byxxx = zeros(maxMat,nx,nx,nx);
E1_byssx = zeros(maxMat,nx);
for i=1:maxMat
    E1_by0(i,:) = model.by0(i,1);
    [E1_byx(i,:),E1_byxx(i,:,:),E1_byss(i,:),E1_byxxx(i,:,:,:),E1_byssx(i,:)] ...
        = CondMoments_Mom3_log(model.byx,model.byxx,model.byss,...
        model.byxxx,model.byssx,model.bysss,...
        model.hx,model.hxx,model.hss,...
        model.hxxx,model.hssx,model.hsss,model.eta,...
        vectorMom3,i,1,setupFilter.appOrder);
end

% Decomposing news for the yield curve states in xSimusing simulated states
%xSim = outEstim.outFilter.outSim.xAll_cu;
T                = size(xSim,2);
nx1              = model.nx1;
E_infl_maxMat_cu = zeros(maxMat+1,T);
E_r_maxMat_cu    = zeros(maxMat+1,T);
yields_cu        = zeros(maxMat,T);
E1_yields_cu      = zeros(maxMat,T);
for t = 1:T
    if setupFilter.pruningOn == 1
        xf_cu = xSim(1:nx,t);
        if setupFilter.appOrder == 1
            E_infl_maxMat_cu(:,t) = E_pai0 + E_paix*xf_cu;
            E_r_maxMat_cu(:,t)    = E_r0   + E_rx*xf_cu;
            yields_cu(:,t)        = model.by0 + model.byx*xf_cu;
            E1_yields_cu(:,t)     = E1_by0 + E1_byx*xf_cu;
        elseif setupFilter.appOrder == 2
            xs_cu = zeros(nx,1);
            xs_cu(1:nx1,1) = xSim(nx+1:end,t);
            xfxf_cu        = kron(xf_cu,xf_cu);
            E_infl_maxMat_cu(:,t) = E_pai0 + E_paix*(xf_cu+xs_cu) + 0.5*reshape(E_paixx,maxMat+1,nx^2)*xfxf_cu + 0.5*E_paiss;
            E_r_maxMat_cu(:,t)    = E_r0   + E_rx*(xf_cu+xs_cu)   + 0.5*reshape(E_rxx,maxMat+1,nx^2)*xfxf_cu + 0.5*E_rss;
            yields_cu(:,t)        = model.by0 + model.byx*(xf_cu+xs_cu) + model.Byxx*xfxf_cu + 0.5*model.byss;
            E1_yields_cu(:,t)     = E1_by0 + E1_byx*(xf_cu+xs_cu) + 0.5*reshape(E1_byxx,maxMat,nx^2)*xfxf_cu + 0.5*E1_byss;
        elseif setupFilter.appOrder == 3
            xs_cu = zeros(nx,1);
            xs_cu(1:nx1,1) = xSim(nx+1:nx+nx1,t);
            xrd_cu = zeros(nx,1);
            xrd_cu(1:nx1,1) = xSim(nx+nx1+1:end,t);
            xfxf_cu         = kron(xf_cu,xf_cu);
            xfxfxf_cu       = kron(xfxf_cu,xf_cu);
            xfxs_cu         = kron(xf_cu,xs_cu);
            E_infl_maxMat_cu(:,t)= E_pai0 + E_paix*(xf_cu+xs_cu + xrd_cu) + 0.5*reshape(E_paixx,maxMat+1,nx^2)*(xfxf_cu + 2*xfxs_cu) ...
                + 1/6*reshape(E_paixxx,maxMat+1,nx^3)*xfxfxf_cu + 0.5*E_paiss + 3/6*E_paissx*xf_cu;
            E_r_maxMat_cu(:,t)    = E_r0   + E_rx*(xf_cu+xs_cu + xrd_cu) + 0.5*reshape(E_rxx,maxMat+1,nx^2)*(xfxf_cu + 2*xfxs_cu) ...
                + 1/6*reshape(E_rxxx,maxMat+1,nx^3)*xfxfxf_cu + 0.5*E_rss + 3/6*E_rssx*xf_cu;
            yields_cu(:,t)        = model.by0 + model.byx*(xf_cu+xs_cu + xrd_cu) + model.Byxx*(xfxf_cu + 2*xfxs_cu) ...
                + model.Byxxx*xfxfxf_cu  + 0.5*model.byss + 3/6*model.byssx*xf_cu;
            E1_yields_cu(:,t)     = E1_by0 + E1_byx*(xf_cu+xs_cu + xrd_cu) + 0.5*reshape(E1_byxx,maxMat,nx^2)*(xfxf_cu + 2*xfxs_cu) ...
                + 1/6*reshape(E1_byxxx,maxMat,nx^3)*xfxfxf_cu + 0.5*E1_byss + 3/6*E1_byssx*xf_cu;
        end
    else
        x_cu   = xSim(1:nx,t);
        xx_cu  = kron(x_cu,x_cu);
        xxx_cu = kron(xx_cu,x_cu);
        E_infl_maxMat_cu(:,t)= E_pai0 + E_paix*x_cu + reshape(E_paixx,maxMat+1,nx^2)*xx_cu + reshape(E_paixxx,maxMat+1,nx^3)*xxx_cu + 0.5*E_paiss + 3/6*E_paissx*x_cu;
        E_r_maxMat_cu(:,t)   = E_r0   + E_rx*x_cu   + reshape(E_rxx,maxMat+1,nx^2)*xx_cu   + reshape(E_rxxx,maxMat+1,nx^3)*xxx_cu + 0.5*E_rss + 3/6*E_rssx*x_cu;
        yields_cu(:,t)       = model.by0 + model.byx*x_cu + model.Byxx*xx_cu + model.Byxxx*xxx_cu + 0.5*model.byss + 3/6*model.byssx*x_cu;
        E1_yields_cu(:,t)    = E1_by0 + E1_byx*x_cu + reshape(E1_byxx,maxMat,nx^2)*xx_cu + reshape(E1_byxxx,maxMat,nx^3)*xxx_cu + 0.5*E1_byss + 3/6*E1_byssx*x_cu;
    end
end
E_realRates_cu      = E_r_maxMat_cu - E_infl_maxMat_cu;

%                     Et[infl_t+m] -                  - Et-1[infl_t+m+1]  
inflNews_per_period = (E_infl_maxMat_cu(1:maxMat,2:T) - E_infl_maxMat_cu(2:maxMat+1,1:T-1))*4;
inflNews            = repmat((1./(1:maxMat)'),1,T-1).*cumsum(inflNews_per_period,1);

realNews_per_period = (E_realRates_cu(1:maxMat,2:T) - E_realRates_cu(2:maxMat+1,1:T-1))*4;
realNews            = repmat((1./(1:maxMat)'),1,T-1).*cumsum(realNews_per_period,1);

yieldsNews          = (yields_cu(:,2:T)- E1_yields_cu(:,1:end-1))*4;
%yieldsNews         = (yields_cu(:,2:T)- yields_cu(:,1:T-1))*4;    %The martingale specification 
% Adding measurement errors to bond yields news
yieldsNews          = yieldsNews+model.params.stdMr20*randn(size(yieldsNews));

premNews            = yieldsNews - inflNews - realNews;

% Output
if T < 1000
    out.realNews           = realNews;
    out.inflNews           = inflNews;
    out.premNews           = premNews;
    out.yieldsNews         = yieldsNews;
end

out.stdRealNewsPct     = std(realNews,0,2)*100;
out.stdInflNewsPct     = std(inflNews,0,2)*100;
out.stdPremNewsPct     = std(premNews,0,2)*100;
out.stdYieldsNewsPct   = std(yieldsNews,0,2)*100;


out.covRealInfl        = zeros(maxMat,1);
out.covRealPrem        = zeros(maxMat,1);
out.covInflPrem        = zeros(maxMat,1);
out.varReal            = out.stdRealNewsPct.^2;
out.varInfl            = out.stdInflNewsPct.^2;
out.varPerm            = out.stdPremNewsPct.^2;
out.varYield           = out.stdYieldsNewsPct.^2;
for i=1:maxMat
    tmp = cov([realNews(i,:)',inflNews(i,:)']);
    out.covRealInfl(i,1) = 2*tmp(1,2)*100^2;
    tmp = cov([realNews(i,:)',premNews(i,:)']);
    out.covRealPrem(i,1) = 2*tmp(1,2)*100^2;
    tmp = cov([inflNews(i,:)',premNews(i,:)']);
    out.covInflPrem(i,1) = 2*tmp(1,2)*100^2;
end
%out.varYieldCheck  = out.varReal + out.varInfl + out.varPerm + out.covRealInfl + out.covRealPrem + out.covInflPrem - out.varYield

out.RealVarRatio       = out.stdRealNewsPct.^2./out.stdYieldsNewsPct.^2;
out.InflVarRatio       = out.stdInflNewsPct.^2./out.stdYieldsNewsPct.^2;
out.PremVarRatio       = out.stdPremNewsPct.^2./out.stdYieldsNewsPct.^2;

matSelect    = [1 3 5 10]*4;
dataMoments  = [0.22 0.21 0.20 nan nan 0.46 0.44 0.46 nan nan 0.23 0.22 0.19 nan];
modelMoments = [ out.stdInflNewsPct(matSelect,1)' nan out.stdYieldsNewsPct(matSelect,1)' nan out.InflVarRatio(matSelect,1)'];
out.summaryTable   = [dataMoments' modelMoments'];

end

